Introduction

Reads were first de-multiplexed with ea-utils (fastq-multx) with an allowed barcode mismatch of 1 (Default). De-multiplexed reads were then subjected to an alignment against the UniVec vector contamination database using the adapter removal tool within ea-utils (fastq-mcf). Barcodes were extracted by remove the last 12 nucleotides from each read.

Next reads were subjected to an alignment against the Oligo library using a DNA aligner. We used Bowtie2…

Data import & preprocessing

We are going to first import the merged counts tables that is outputted from the Nextflow pipeline. This table already has the raw counts summarised from the alignments and merged into a single matrix.

Import data

# Import CD4 counts
counts = readRDS('../data/CD4-counts.rda')
head(counts)
##         cd4_bcell_alone_1 cd4_cd8_mansc1_1 cd4_cd8_mansc1_2 cd4_cd8_mansc1_3
## oligo_1              9153               45               53               38
## oligo_2             15918               76               94               80
## oligo_3              7795               28               29               34
## oligo_4              9068              132               39               44
## oligo_5              8207               25               55               41
## oligo_6              5805               23               42               32
##         cd4_cd8_mock_1 cd4_cd8_mock_2 cd4_cd8_snord73_1 cd4_cd8_snord73_2
## oligo_1             39             55                25                47
## oligo_2             88             88                35               108
## oligo_3             25             26                14                42
## oligo_4             46             36                23                43
## oligo_5             45             42                19                54
## oligo_6             35             35                11                37
##         cd4_library_maxi_1 cd4_library_maxi_2 cd4_mansc1_1 cd4_mansc1_2
## oligo_1               1969               1955         5790         5149
## oligo_2               3357               2816        10287         8929
## oligo_3               1506               1595         4868         4243
## oligo_4               1604               1761         4565         3622
## oligo_5               1226               1334         4701         3600
## oligo_6               1251               1191         3228         3239
##         cd4_mansc1_3 cd4_mock_1 cd4_mock_2 cd4_snord73_1 cd4_snord73_2
## oligo_1         5117       6086       1844          6762          7154
## oligo_2         9485      10287       2803         13097         12108
## oligo_3         4508       5181       1422          6088          5464
## oligo_4         4399       5286       1396          5650          6561
## oligo_5         3825       4086       1275          5025          4895
## oligo_6         3035       3849       1305          4501          3700
##         cd4_snord73_3
## oligo_1          7346
## oligo_2         13418
## oligo_3          5936
## oligo_4          6067
## oligo_5          5383
## oligo_6          4226
# Import CD4 metadata
metadata = readRDS('../data/CD4-metadata.rda')
head(metadata)
##                                Sample Group Barcode
## cd4_library_maxi_1 cd4_library_maxi_1   DNA GGAATTC
## cd4_library_maxi_2 cd4_library_maxi_2   DNA TTAGTGG
## cd4_bcell_alone_1   cd4_bcell_alone_1 Bcell ACAACGA
## cd4_mock_1                 cd4_mock_1  Mock ACCAATG
## cd4_mock_2                 cd4_mock_2  Mock CAACAGC
## cd4_snord73_1           cd4_snord73_1 SNORD CACACGT
# Import MultiQC results
report = readRDS('../data/CD4-multiqc.rda')
head(report)
## # A tibble: 6 × 47
##   Sample   raw_total_seque… filtered_sequen… sequences is_sorted `1st_fragments`
##   <chr>               <dbl>            <dbl>     <dbl>     <dbl>           <dbl>
## 1 cd4_bce…         32649199                0  32649199         1        32649199
## 2 cd4_cd8…           204206                0    204206         1          204206
## 3 cd4_cd8…           196064                0    196064         1          196064
## 4 cd4_cd8…           200980                0    200980         1          200980
## 5 cd4_cd8…           209112                0    209112         1          209112
## 6 cd4_cd8…           226005                0    226005         1          226005
## # … with 41 more variables: last_fragments <dbl>, reads_mapped <dbl>,
## #   reads_mapped_and_paired <dbl>, reads_unmapped <dbl>,
## #   reads_properly_paired <dbl>, reads_paired <dbl>, reads_duplicated <dbl>,
## #   reads_MQ0 <dbl>, reads_QC_failed <dbl>, `non-primary_alignments` <dbl>,
## #   total_length <dbl>, total_first_fragment_length <dbl>,
## #   total_last_fragment_length <dbl>, bases_mapped <dbl>,
## #   `bases_mapped_(cigar)` <dbl>, bases_trimmed <dbl>, …
# Import oligo annotations
annotations = readRDS('../data/CD4-annotations.rda')
head(annotations)
## # A tibble: 6 × 3
##   oligo_id gene_name gene_id  
##   <chr>    <chr>     <chr>    
## 1 oligo_1  MAGEA1    MAGEA1 #1
## 2 oligo_2  <NA>      NA #1    
## 3 oligo_3  <NA>      NA #2    
## 4 oligo_4  <NA>      NA #3    
## 5 oligo_5  <NA>      NA #4    
## 6 oligo_6  <NA>      NA #5

Mapping quality

# Mean
mean(report$raw_total_sequences)
## [1] 11596529
median(report$raw_total_sequences)
## [1] 7196006
range(report$raw_total_sequences)
## [1]   114571 32649199

Number of mapped reads

p1 = report %>% 
  select(Sample, reads_unmapped, reads_mapped) %>% 
  mutate(Sample = str_replace(Sample, '-bamstats', '')) %>% 
  pivot_longer(cols = c("reads_unmapped", "reads_mapped"), names_to = "feature", values_to = "reads") %>% 
    mutate(feature = if_else(feature == "reads_unmapped", "Unmapped", feature)) %>% 
  mutate(feature = if_else(feature == "reads_mapped", "Mapped", feature)) %>% 
  ggplot(aes(x = reorder(Sample, reads), y = log10(reads), fill = feature)) +
  geom_col() +
  coord_flip() +
  theme_classic(base_size = 10) +
  theme(legend.position = "top") +
  scale_fill_brewer(palette = "Set2") +
  xlab("") +
  ylab("Number of reads (log10)") +
  ggtitle("CD4: mapping rate")
p1

### Percent of mapped reads
p2 = report %>% 
  select(Sample, reads_unmapped_percent, reads_mapped_percent ) %>% 
  mutate(Sample = str_replace(Sample, '-bamstats', '')) %>% 
  pivot_longer(cols = c("reads_unmapped_percent", "reads_mapped_percent"), names_to = "feature", values_to = "reads") %>% 
  mutate(reads = reads / 100) %>% 
  mutate(feature = if_else(feature == "reads_unmapped_percent", "Unmapped", feature)) %>% 
  mutate(feature = if_else(feature == "reads_mapped_percent", "Mapped", feature)) %>% 
  ggplot(aes(x = reorder(Sample, reads), y = reads, fill = feature)) +
  geom_col() +
  coord_flip() +
  theme_classic(base_size = 10) +
  theme(legend.position = "top") +
  scale_fill_brewer(palette = "Set2") +
  scale_y_continuous(labels = scales::percent) +
  xlab("") +
  ylab("Percent of reads")
p2

# Plot together
p = (p1 | p2) + plot_layout(guides = "collect") & theme(legend.position = "bottom")
p

ggsave("figures/Benchmark-CD4-mapping-qc.pdf", p, width = 8, height = 5)

Normalize data

# Match the sample ID's between samples
counts = counts[,match(metadata$Sample, colnames(counts))]

# Make DESeq2 object
ddsMat <- DESeqDataSetFromMatrix(countData = counts,
                                 colData = metadata,
                                 design = ~Group)

# Get normalize counts
ddsMat <- estimateSizeFactors(ddsMat)

# Filter flow counts
ddsMat.fil = ddsMat %>% 
  filter_low(percentile = 0.05)

# Get normalized counts 
counts.norm.fil <- (counts(ddsMat.fil, normalized = TRUE) + 0.5) %>% 
  as.data.frame() %>% 
  rownames_to_column("oligo_id")

# Create a long table 
counts.norm.fil.long = counts.norm.fil %>% 
  pivot_longer(cols = -oligo_id, names_to = "sampleId", values_to = "norm")

Mock replicates

counts.norm.fil %>% 
  ggplot(aes(x = log10(cd4_mock_1), y = log10(cd4_mock_2))) +
  geom_point(colour = alpha("grey", 0.7)) +
  xlab("Mock replicate #1 (log10)") +
  ylab("Mock replicate #2 (log10)") +
  theme_light() +
  stat_cor(method = "pearson", aes(label = ..r.label..)) +
  geom_smooth(method = "lm", se = FALSE)
## `geom_smooth()` using formula 'y ~ x'

### SNORD73 replicates

counts.norm.fil %>% 
  ggplot(aes(x = log10(cd4_snord73_1), y = log10(cd4_snord73_2))) +
  geom_point(colour = alpha("grey", 0.7)) +
  xlab("SNORD73 replicate #1 (log10)") +
  ylab("SNORD73 replicate #2 (log10)") +
  theme_light(base_size = 11) +
  stat_cor(method = "pearson", aes(label = ..r.label..)) +
  geom_smooth(method = "lm", se = FALSE)
## `geom_smooth()` using formula 'y ~ x'

counts.norm.fil %>% 
  ggplot(aes(x = log10(cd4_snord73_2), y = log10(cd4_snord73_3))) +
  geom_point(colour = alpha("grey", 0.7)) +
  xlab("SNORD73 replicate #1 (log10)") +
  ylab("SNORD73 replicate #2 (log10)") +
  theme_light(base_size = 11) +
  stat_cor(method = "pearson", aes(label = ..r.label..)) +
  geom_smooth(method = "lm", se = FALSE) 
## `geom_smooth()` using formula 'y ~ x'

counts.norm.fil %>% 
  ggplot(aes(x = log10(cd4_snord73_1), y = log10(cd4_snord73_2))) +
  geom_point(colour = alpha("grey", 0.7)) +
  xlab("SNORD73 replicate #1 (log10)") +
  ylab("SNORD73 replicate #2 (log10)") +
  theme_light(base_size = 11) +
  stat_cor(method = "pearson", aes(label = ..r.label..)) +
  geom_smooth(method = "lm", se = FALSE)
## `geom_smooth()` using formula 'y ~ x'

MANSC1 replicates

counts.norm.fil %>% 
  ggplot(aes(x = log10(cd4_mansc1_1), y = log10(cd4_mansc1_2))) +
  geom_point(colour = alpha("grey", 0.7)) +
  xlab("MANSC1 replicate #1 (log10)") +
  ylab("MANSC1 replicate #2 (log10)") +
  theme_light(base_size = 11) +
  stat_cor(method = "pearson", aes(label = ..r.label..)) +
  geom_smooth(method = "lm", se = FALSE) 
## `geom_smooth()` using formula 'y ~ x'

counts.norm.fil %>% 
  ggplot(aes(x = log10(cd4_mansc1_2), y = log10(cd4_mansc1_3))) +
  geom_point(colour = alpha("grey", 0.7)) +
  xlab("MANSC1 replicate #2 (log10)") +
  ylab("MANSC1 replicate #3 (log10)") +
  theme_light(base_size = 11) +
  stat_cor(method = "pearson", aes(label = ..r.label..)) +
  geom_smooth(method = "lm", se = FALSE) 
## `geom_smooth()` using formula 'y ~ x'

counts.norm.fil %>% 
  ggplot(aes(x = log10(cd4_mansc1_1), y = log10(cd4_mansc1_3))) +
  geom_point(colour = alpha("grey", 0.7)) +
  xlab("MANSC1 replicate #1 (log10)") +
  ylab("MANSC1 replicate #3 (log10)") +
  theme_light(base_size = 11) +
  stat_cor(method = "pearson", aes(label = ..r.label..)) +
  geom_smooth(method = "lm", se = FALSE) 
## `geom_smooth()` using formula 'y ~ x'

Plot normalized histogram

counts.norm.fil.long %>% 
  ggplot(aes(x = log10(norm))) +
  geom_histogram(aes(y=..density..), binwidth=0.25, colour="black", fill="white") +
  geom_density(alpha = 0.2, fill = "#FF6666") +
  facet_wrap(sampleId~.) +
  theme_light(base_size = 11) +
  ylab("Density") +
  xlab("Normalized reads per oligo (log10)") 


Determine dropout oligo’s

SNORD73 vs. Mock

snord_mock.res = compute_stats(ddsMat.fil, trt = "SNORD", ref = "Mock", fitType = "parametric", lfcThreshold = 0.25, alpha = 0.05) %>% 
  left_join(annotations) %>% 
  mutate(log2FoldChange = round(log2FoldChange, 2)) %>% 
  mutate(log2FoldChangeShrunken = round(log2FoldChangeShrunken, 2)) %>% 
  mutate(gene_name = if_else(gene_name == "RPS3Amut", "SNORDmut", gene_name)) %>% 
  mutate(gene_name = if_else(gene_name == "RPS3Awt", "SNORDwt", gene_name))
## using pre-existing size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## -- note: fitType='parametric', but the dispersion trend was not well captured by the
##    function: y = a/x + b, and a local regression fit was automatically substituted.
##    specify fitType='local' or 'mean' to avoid this message next time.
## final dispersion estimates
## fitting model and testing
## using 'normal' for LFC shrinkage, the Normal prior from Love et al (2014).
## 
## Note that type='apeglm' and type='ashr' have shown to have less bias than type='normal'.
## See ?lfcShrink for more details on shrinkage type, and the DESeq2 vignette.
## Reference: https://doi.org/10.1093/bioinformatics/bty895
## 
## out of 4525 with nonzero total read count
## adjusted p-value < 0.05
## LFC > 0.25 (up)    : 0, 0%
## LFC < -0.25 (down) : 2, 0.044%
## outliers [1]       : 0, 0%
## low counts [2]     : 0, 0%
## (mean count < 14)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results
## Joining, by = "oligo_id"
# Significant oligos
snord_mock.res %>% 
  filter(padj < 0.05)
##     oligo_id baseMean log2FoldChange     lfcSE       stat       pvalue
## 1 oligo_1889 195.6394          -2.50 0.2088072 -10.762587 2.584900e-27
## 2 oligo_4271 162.4496          -2.35 0.2284288  -9.182969 2.096831e-20
##           padj log2FoldChangeShrunken gene_name     gene_id
## 1 1.169667e-23                  -0.83  SNORDmut RPS3Amut #1
## 2 4.744080e-17                  -0.69  SNORDmut RPS3Amut #2
# Browse
snord_mock.res %>% 
  select(gene_id, log2FoldChange, padj, gene_name, gene_id) %>% 
  DT::datatable(extensions = 'Buttons', options = list(
    dom = 'Bfrtip',
    buttons = c('copy', 'csv', 'excel', 'pdf', 'print')
  ))
# Plot MA-plot
snord_mock.maplot = snord_mock.res %>% 
  ggplot(aes(x = log2(baseMean), y = log2FoldChange, label = gene_name)) +
  geom_point(colour = alpha("grey", 0.7)) +
  geom_point(colour = alpha("grey", 0.7)) +
  geom_point(data = filter(snord_mock.res, gene_name %in% c("SNORDmut")), size = 3, colour = "#efc383") +
  geom_point(data = filter(snord_mock.res, gene_name %in% c("SNORDwt")), size = 3, colour = "#a6a6a6") +
  geom_text_repel(data = filter(snord_mock.res, gene_name %in% c("SNORDmut", "SNORDwt")), box.padding = 0.5, max.overlaps = Inf, size = 2.5) +
  theme_classic() +
  theme_classic() +
  xlab("Mean counts (log2)") +
  ylab("Fold change log2(SNORD73/Mock)")  +
  ggtitle("SNORD")
ggplotly(snord_mock.maplot)

MANSC1 vs. Mock

mansc_mock.res = compute_stats(ddsMat.fil, trt = "MANSC", ref = "Mock", fitType = "parametric", lfcThreshold = 0.25, alpha = 0.05) %>% 
  left_join(annotations) %>% 
  mutate(log2FoldChange = round(log2FoldChange, 2)) %>% 
  mutate(log2FoldChangeShrunken = round(log2FoldChangeShrunken, 2)) %>% 
  mutate(gene_name = if_else(gene_name == "RPS3Amut", "SNORDmut", gene_name)) %>% 
  mutate(gene_name = if_else(gene_name == "RPS3Awt", "SNORDwt", gene_name))
## using pre-existing size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## -- note: fitType='parametric', but the dispersion trend was not well captured by the
##    function: y = a/x + b, and a local regression fit was automatically substituted.
##    specify fitType='local' or 'mean' to avoid this message next time.
## final dispersion estimates
## fitting model and testing
## using 'normal' for LFC shrinkage, the Normal prior from Love et al (2014).
## 
## Note that type='apeglm' and type='ashr' have shown to have less bias than type='normal'.
## See ?lfcShrink for more details on shrinkage type, and the DESeq2 vignette.
## Reference: https://doi.org/10.1093/bioinformatics/bty895
## 
## out of 4525 with nonzero total read count
## adjusted p-value < 0.05
## LFC > 0.25 (up)    : 0, 0%
## LFC < -0.25 (down) : 2, 0.044%
## outliers [1]       : 0, 0%
## low counts [2]     : 0, 0%
## (mean count < 10)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results
## Joining, by = "oligo_id"
# Significant oligos
mansc_mock.res %>% 
  filter(padj < 0.05) 
##     oligo_id baseMean log2FoldChange     lfcSE       stat       pvalue
## 1 oligo_1895 618.9236          -1.94 0.1207156 -13.993815 8.502020e-45
## 2 oligo_4277 429.8109          -2.05 0.3513332  -5.119629 1.530689e-07
##           padj log2FoldChangeShrunken gene_name      gene_id
## 1 3.847164e-41                  -1.21 MANSC1mut MANSC1mut #1
## 2 3.463184e-04                  -0.33 MANSC1mut MANSC1mut #2
# Browse
mansc_mock.res %>% 
  select(gene_id, log2FoldChange, padj, gene_name, gene_id) %>% 
  DT::datatable(extensions = 'Buttons', options = list(
    dom = 'Bfrtip',
    buttons = c('copy', 'csv', 'excel', 'pdf', 'print')
  ))
# Plot MA-plot
mansc_mock.maplot = mansc_mock.res %>% 
  ggplot(aes(x = log2(baseMean), y = log2FoldChange, label = gene_name)) +
  geom_point(colour = alpha("grey", 0.7)) +
  geom_point(colour = alpha("grey", 0.7)) +
  geom_point(data = filter(mansc_mock.res, gene_name %in% c("MANSC1mut")), size = 3, colour = "#aec3a0") +
  geom_point(data = filter(mansc_mock.res, gene_name %in% c("MANSC1wt")), size = 3, colour = "#a6a6a6") +
  geom_text_repel(data = filter(mansc_mock.res, gene_name %in% c("MANSC1mut", "MANSC1wt")), box.padding = 0.5, max.overlaps = Inf, size = 2.5) +
  theme_classic() +
  theme_classic() +
  xlab("Mean counts (log2)") +
  ylab("Fold change log2(MANSC1/Mock)")  +
  ggtitle("MANSC1")
ggplotly(mansc_mock.maplot)

Plot together

p = (snord_mock.maplot + ggtitle("") | mansc_mock.maplot + ggtitle("")) & ylim(-3, 0.75) & xlim(7, 11)
p

ggsave("figures/CD4-benchmark.pdf", p, width = 6, height = 4)

MANSC1 vs. SNORD73

mansc_snord.res = compute_stats(ddsMat.fil, trt = "MANSC", ref = "SNORD", fitType = "parametric", lfcThreshold = 0.25, alpha = 0.05) %>% 
  left_join(annotations) %>% 
  mutate(log2FoldChange = round(log2FoldChange, 2)) %>% 
  mutate(log2FoldChangeShrunken = round(log2FoldChangeShrunken, 2)) %>% 
  mutate(gene_name = if_else(gene_name == "RPS3Amut", "SNORDmut", gene_name)) %>% 
  mutate(gene_name = if_else(gene_name == "RPS3Awt", "SNORDwt", gene_name))

# Get significant hits
mansc_snord.res %>%  
  filter(padj < 0.05) 
mansc1_snord.res

# Browse
mansc_snord.res %>% 
  select(gene_id, log2FoldChange, padj, gene_name, gene_id) %>% 
  DT::datatable(extensions = 'Buttons', options = list(
    dom = 'Bfrtip',
    buttons = c('copy', 'csv', 'excel', 'pdf', 'print')
  ))

# Plot MA-plot
mansc_snord.maplot = mansc_snord.res %>% 
  ggplot(aes(x = log2(baseMean), y = log2FoldChange, label = gene_name)) +
  geom_point(colour = alpha("grey", 0.7)) +
  geom_point(colour = alpha("grey", 0.7)) +
  geom_point(data = filter(mansc_snord.res, gene_name %in% c("MANSC1mut")), size = 3, colour = "#aec3a0") +
  geom_point(data = filter(mansc_snord.res, gene_name %in% c("MANSC1wt")), size = 3, colour = "#a6a6a6") +
  geom_text_repel(data = filter(mansc_snord.res, gene_name %in% c("MANSC1mut", "MANSC1wt")), box.padding = 0.5, max.overlaps = Inf, size = 2.5) +
  theme_classic() +
  theme_classic() +
  xlab("Mean counts (log2)") +
  ylab("Fold change log2(MANSC1/Mock)") 
ggplotly(mansc_snord.maplot)

Session info

sessionInfo()
## R version 4.1.3 (2022-03-10)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.5 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/liblapack.so.3
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] stats4    stats     graphics  grDevices datasets  utils     methods  
## [8] base     
## 
## other attached packages:
##  [1] forcats_0.5.1               stringr_1.4.0              
##  [3] dplyr_1.0.9                 purrr_0.3.4                
##  [5] readr_2.1.2                 tidyr_1.2.0                
##  [7] tibble_3.1.7                tidyverse_1.3.1            
##  [9] ggpubr_0.4.0                plotly_4.10.0              
## [11] ggplotify_0.1.0             readxl_1.4.0               
## [13] DESeq2_1.34.0               SummarizedExperiment_1.24.0
## [15] Biobase_2.54.0              MatrixGenerics_1.6.0       
## [17] matrixStats_0.62.0          GenomicRanges_1.46.1       
## [19] GenomeInfoDb_1.30.1         IRanges_2.28.0             
## [21] S4Vectors_0.32.4            BiocGenerics_0.40.0        
## [23] ggrepel_0.9.1               ggplot2_3.3.6              
## [25] ggsci_2.9                   patchwork_1.1.1            
## [27] knitr_1.39                 
## 
## loaded via a namespace (and not attached):
##   [1] colorspace_2.0-3       ggsignif_0.6.3         ellipsis_0.3.2        
##   [4] XVector_0.34.0         fs_1.5.2               rstudioapi_0.13       
##   [7] farver_2.1.0           DT_0.25                bit64_4.0.5           
##  [10] AnnotationDbi_1.56.2   fansi_1.0.3            lubridate_1.8.0       
##  [13] xml2_1.3.3             splines_4.1.3          cachem_1.0.6          
##  [16] geneplotter_1.72.0     jsonlite_1.8.0         broom_0.8.0           
##  [19] annotate_1.72.0        dbplyr_2.1.1           png_0.1-7             
##  [22] BiocManager_1.30.18    compiler_4.1.3         httr_1.4.2            
##  [25] backports_1.4.1        assertthat_0.2.1       Matrix_1.4-0          
##  [28] fastmap_1.1.0          lazyeval_0.2.2         cli_3.3.0             
##  [31] htmltools_0.5.2        tools_4.1.3            gtable_0.3.0          
##  [34] glue_1.6.2             GenomeInfoDbData_1.2.7 Rcpp_1.0.8.3          
##  [37] carData_3.0-5          cellranger_1.1.0       jquerylib_0.1.4       
##  [40] vctrs_0.4.1            Biostrings_2.62.0      nlme_3.1-155          
##  [43] crosstalk_1.2.0        xfun_0.30              rvest_1.0.2           
##  [46] lifecycle_1.0.1        renv_0.15.4            rstatix_0.7.0         
##  [49] XML_3.99-0.10          zlibbioc_1.40.0        scales_1.2.0          
##  [52] hms_1.1.1              parallel_4.1.3         RColorBrewer_1.1-3    
##  [55] yaml_2.3.5             memoise_2.0.1          yulab.utils_0.0.5     
##  [58] sass_0.4.1             stringi_1.7.6          RSQLite_2.2.15        
##  [61] highr_0.9              genefilter_1.76.0      BiocParallel_1.28.3   
##  [64] rlang_1.0.2            pkgconfig_2.0.3        bitops_1.0-7          
##  [67] evaluate_0.15          lattice_0.20-45        labeling_0.4.2        
##  [70] htmlwidgets_1.5.4      bit_4.0.4              tidyselect_1.1.2      
##  [73] magrittr_2.0.3         R6_2.5.1               generics_0.1.2        
##  [76] DelayedArray_0.20.0    DBI_1.1.2              mgcv_1.8-39           
##  [79] pillar_1.7.0           haven_2.5.0            withr_2.5.0           
##  [82] survival_3.2-13        KEGGREST_1.34.0        abind_1.4-5           
##  [85] RCurl_1.98-1.6         modelr_0.1.8           crayon_1.5.1          
##  [88] car_3.0-13             utf8_1.2.2             tzdb_0.3.0            
##  [91] rmarkdown_2.14         locfit_1.5-9.6         grid_4.1.3            
##  [94] data.table_1.14.2      blob_1.2.3             reprex_2.0.1          
##  [97] digest_0.6.29          xtable_1.8-4           gridGraphics_0.5-1    
## [100] munsell_0.5.0          viridisLite_0.4.0      bslib_0.3.1